Méthodes de Partionnement et d'apprentissage non supervisé

Classification Hiérarchique et Kmeans

Anne Badel et Frédéric Guyon

2019-02-19

Partionnement et apprentissage

Partionnement=Clustering

Y a-t-il des groupes ?

Y a-t-il des groupes ?

Partionnement=Clustering

Oui, 4 groupes.

Oui, 4 groupes.

Apprentissage

2 groupes

Apprentissage: Séparation linéaire

2 groupes

Méthodes

Trois grands principes de méthodes basées sur:

En fait, trois façons de voir les mêmes algorithmes

Géométrie et distances

On considère les données comme des points de \(R^n\):

Les données

Ces données sont un classique des méthodes d'apprentissage

Dans un premier temps, regardons les données

dim(mes.iris)
[1] 150   4
head(mes.iris)
  Sepal.Length Sepal.Width Petal.Length Petal.Width
1          5.1         3.5          1.4         0.2
2          4.9         3.0          1.4         0.2
3          4.7         3.2          1.3         0.2
4          4.6         3.1          1.5         0.2
5          5.0         3.6          1.4         0.2
6          5.4         3.9          1.7         0.4

Visualisation des données

On peut ensuite essayer de visualiser les données

plot(mes.iris)

Géométrie et distances

Sur la base d'une distance (souvent euclidienne)

Distances

Définition d'une distance : fonction positive de deux variables

  1. d(x,y) \(\ge\) 0
  2. d(x,y)=d(y,x)
  3. d(x,y)=0 \(\Longleftrightarrow\) x=y
  4. inégalité triangulaire: d(x,z) \(\le\) d(x,y)+d(y,z)

Si 1,2,3 : dissimilarité

Distances utilisées dans R

Distances

Distances utilisées dans R

Autres distances non géométriques (pour information)

Utilisées en bio-informatique:

\[d("BONJOUR", "BONSOIR")=2\]

Distances plus classiques en génomiques

Comme vu lors de la séance 3, il y a d'autres mesures de distances :

ne sont pas des distances, mais indices de dissimilarité:

rq : lors du TP, sur les données d'expression RNA-seq, nous utiliserons le coefficient de corrélation de Spearman et la distance associée, \(d = 1-r^2\)

Distances entre groupes

\(D(C_1,C_2) = \min_{i \in C_1, j \in C_2} D(x_i, x_j)\)

\(D(C_1,C_2) = \max_{i \in C_1, j \in C_2} D(x_i, x_j)\)

\(D(C_1,C_2) = \frac{1}{N_1 N_2} \sum_{i \in C_1, j \in C_2} D(x_i, x_j)\)

\(d^2(C_i,C_j)= I_{intra}(C_i \cup C_j)-I_{intra}(C_i)-I_{intra}(C_j)\) \(D(C_1,C_2) = \sqrt{\frac{N_1N_2}{N_1 + N_2}} \| m_1 -m_2 \|\)

Distances entre groupes

Distances entre groupes

Les données

Ces données sont un classique des méthodes d'apprentissage

Dans un premier temps, regardons les données

dim(mes.iris)
[1] 150   4
head(mes.iris)
  Sepal.Length Sepal.Width Petal.Length Petal.Width
1          5.1         3.5          1.4         0.2
2          4.9         3.0          1.4         0.2
3          4.7         3.2          1.3         0.2
4          4.6         3.1          1.5         0.2
5          5.0         3.6          1.4         0.2
6          5.4         3.9          1.7         0.4
str(mes.iris)
'data.frame':   150 obs. of  4 variables:
 $ Sepal.Length: num  5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
 $ Sepal.Width : num  3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
 $ Petal.Length: num  1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
 $ Petal.Width : num  0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
summary(mes.iris)
  Sepal.Length    Sepal.Width     Petal.Length    Petal.Width   
 Min.   :4.300   Min.   :2.000   Min.   :1.000   Min.   :0.100  
 1st Qu.:5.100   1st Qu.:2.800   1st Qu.:1.600   1st Qu.:0.300  
 Median :5.800   Median :3.000   Median :4.350   Median :1.300  
 Mean   :5.843   Mean   :3.057   Mean   :3.758   Mean   :1.199  
 3rd Qu.:6.400   3rd Qu.:3.300   3rd Qu.:5.100   3rd Qu.:1.800  
 Max.   :7.900   Max.   :4.400   Max.   :6.900   Max.   :2.500  

Visualisation des données

On peut ensuite essayer de visualiser les données

plot(mes.iris)

Visualisation des données

image(1:nb.var, 1:nb.iris ,t(as.matrix(mes.iris)))

Nettoyage des données (1)

Avant de commencer à travailler, il est nécessaire de commencer par vérifier que :

sum(is.na(mes.iris))
[1] 0

Nettoyage des données (2)

iris.var <- apply(mes.iris, 2, var)
kable(iris.var, digits = 3)
x
Sepal.Length 0.686
Sepal.Width 0.190
Petal.Length 3.116
Petal.Width 0.581
sum(apply(mes.iris, 2, var) == 0)
[1] 0

Normalisation

Afin de pouvoir considérer que toutes les variables sont à la même échelle, il est parfois nécessaire de normaliser les données.

mes.iris.centre <- scale(mes.iris, center=TRUE, scale=FALSE)
mes.iris.scaled <- scale(mes.iris, center=TRUE, scale=TRUE)

On peut visuellement regarder l'effet de la normalisation :

par(mfrow=c(1,2))
plot(mes.iris)

plot(mes.iris.scaled)
par(mfrow=c(1,1))

! ne pas faire si "grosses" données

par(mfrow=c(1,2))
image(1:nb.var, 1:nb.iris, t(as.matrix(mes.iris)), main="données brutes")
image(1:nb.var, 1:nb.iris, t(as.matrix(mes.iris.scaled)), main="données normalisées")

par(mfrow=c(1,1))

par(mfrow=c(1,2))
biplot(prcomp(mes.iris), main="données non normalisées")
biplot(prcomp(mes.iris, scale=T), main="données normalisées")

par(mfrow=c(1,1))

La matrice de distance

Nous utilisons la distance euclidienne

iris.euc <- dist(mes.iris)
iris.scale.euc <- dist(mes.iris.scaled)
par(mfrow=c(1,2))
image(t(as.matrix(iris.euc)), main="données brutes")
image(t(as.matrix(iris.scale.euc)), main="données normalisées")

par(mfrow=c(1,1))

La classification hiérarchique

Principe

Notion importante, cf distances

L'algorithme

étape 1 :

au départ

identification des individus les plus proches

construction du dendrogramme

étape j :

calcul des nouveaux représentants BE et CD

calcul des distances de l'individu restant A aux points moyens

A est plus proche de ...

dendrogramme

pour finir

dendrogramme final

Je ne fais pas attention à ce que je fais

iris.hclust <- hclust(iris.euc)
plot(iris.hclust, hang=-1, cex=0.5)

c'est à dire aux options des fonctions dist et hclust

Sur données normalisées

iris.scale.hclust <- hclust(iris.scale.euc)
plot(iris.scale.hclust, hang=-1, cex=0.5)

par(mfrow=c(1,2))
plot(iris.hclust, hang=-1, cex=0.5)
plot(iris.scale.hclust, hang=-1, cex=0.5)

par(mfrow=c(1,1))

En utilisant une autre métrique

iris.scale.max <- dist(mes.iris.scaled, method="max")
iris.scale.hclust.max <- hclust(iris.scale.max)
par(mfrow=c(1,2))
plot(iris.scale.hclust, hang=-1, cex=0.5)
plot(iris.scale.hclust.max, hang=-1, cex=0.5)

par(mfrow=c(1,1))

En utilisant un autre critère d'aggrégation

iris.scale.hclust.ward <- hclust(iris.scale.euc, method="ward.D2")
par(mfrow=c(1,2))
plot(iris.scale.hclust, hang=-1, cex=0.5)
plot(iris.scale.hclust.ward, hang=-1, cex=0.5)

par(mfrow=c(1,1))

Les k-means

Les individus dans le plan

L'algorithme

étape 1 :

choix des centres provisoires

calcul des distances aux centres provisoires

et affectation à un cluster

calcul des nouveaux centres de classes

étape j :

fin :

arrêt :

Un premier k-means en 5 groupes

iris.scale.kmeans5 <- kmeans(mes.iris.scaled, center=5)
iris.scale.kmeans5
K-means clustering with 5 clusters of sizes 12, 50, 37, 24, 27

Cluster means:
  Sepal.Length Sepal.Width Petal.Length Petal.Width
1    1.6316667  0.06766667    2.5420000  0.85066667
2   -0.8373333  0.37066667   -2.2960000 -0.95333333
3    0.3863964 -0.20598198    1.0095676  0.37363964
4    0.6858333  0.00100000    1.7503333  0.96316667
5   -0.3137037 -0.43511111    0.1827407  0.01918519

Clustering vector:
  [1] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 5 3 5 3 5 3 5 5 5 5 3 5 3 3 5 3 5 3 5 3 3 3 3 3 3 3 5 5 5 5 3 5 3 3 3 5 5 5 3 5 5 5 5 5 3 5 5 4 3 1 4 4 1 5 1 4 1 4 4 4 3 4 4 4 1 1 3 4 3 1 3 4 1 3 3 4 1 1 1 4 3 3 1 4 4 3 4 4 4 3 4 4 4 3
[148] 4 4 3

Within cluster sum of squares by cluster:
[1]  4.655000 15.151000 11.963784  5.462500  9.228889
 (between_SS / total_SS =  93.2 %)

Available components:

[1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss" "betweenss"    "size"         "iter"         "ifault"      
iris.scale.kmeans5$cluster
  [1] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 5 3 5 3 5 3 5 5 5 5 3 5 3 3 5 3 5 3 5 3 3 3 3 3 3 3 5 5 5 5 3 5 3 3 3 5 5 5 3 5 5 5 5 5 3 5 5 4 3 1 4 4 1 5 1 4 1 4 4 4 3 4 4 4 1 1 3 4 3 1 3 4 1 3 3 4 1 1 1 4 3 3 1 4 4 3 4 4 4 3 4 4 4 3
[148] 4 4 3
table(iris.scale.kmeans5$cluster)

 1  2  3  4  5 
12 50 37 24 27 

Visualisation des clusters

plot(iris.scaled.acp, col.ind=iris.scale.kmeans5$cluster, choix="ind")

Combien de clusters ?

Quand une partition est-elle bonne ?

Classification hiérarchique

La coupure de l’arbre à un niveau donné construit une partition. la coupure doit se faire :

plot(iris.scale.hclust.ward, hang=-1, cex=0.5)

K-means

I.intra = numeric(length=10)
I.intra[1] = kmeans(mes.iris.scaled, centers=2)$totss
for (i in 2:10) {
  kmi <- kmeans(mes.iris.scaled, centers=i)
  I.intra[i] <- kmi$tot.withinss
}
plot(1:10, I.intra, type="l")

iris.scale.kmeans3 <- kmeans(mes.iris.scaled, center=3)
plot(iris.scaled.acp, col.ind=iris.scale.kmeans3$cluster, choix="ind")

Heatmap

heatmap(mes.iris.scaled)

my_group=as.numeric(as.factor(substr(variete, 1 , 2)))
my_col=brewer.pal(3, "Set1")[my_group]
heatmap(mes.iris.scaled, RowSideColors=my_col)

Comparaison de clustering: Rand Index

Mesure de similarité entre deux clustering

à partir du nombre de fois que les classifications sont d'accord

\[R=\frac{m+s}{t}\]

Comparaison de clustering: Adjusted Rand Index

\[ ARI=\frac{RI-Expected RI}{Max RI -Expected RI}\]

Comparaison des résultats des deux classifications

cluster.kmeans3 <- iris.scale.kmeans3$cluster
cluster.hclust5 <- cutree(iris.scale.hclust.ward, k=5)
table(cluster.hclust5, cluster.kmeans3)
               cluster.kmeans3
cluster.hclust5  1  2  3
              1 50  0  0
              2  0  2 36
              3  0  0 26
              4  0 24  0
              5  0 12  0
par(mfrow=c(1,2))
plot(iris.scaled.acp, col.ind=cluster.kmeans3, choix="ind", title="kmeans en 3 groupes", cex=0.6)
plot(iris.scaled.acp, col.ind=cluster.hclust5, choix="ind", title="hclust en 5 groupes", cex=0.6)

par(mfrow=c(1,1))

Comparaison avec la réalité

La réalité

variete <- iris[,5]
table(variete)
variete
    setosa versicolor  virginica 
        50         50         50 
plot(iris.scaled.acp, col.ind=variete, choix="ind")

Comparer k-means avec la réalité

conf.kmeans <- table(variete, cluster.kmeans3)
kable(conf.kmeans)
1 2 3
setosa 50 0 0
versicolor 0 2 48
virginica 0 36 14

Setosa vs !Setosa

Visualisation

variete2 <- rep("notSetosa", 150)
variete2[variete=="setosa"] <- "setosa"
variete2 = factor(variete2)
table(variete2)
variete2
notSetosa    setosa 
      100        50 
par(mfrow=c(1,2))
plot(iris.scaled.acp, col.ind=variete2, title="variétés observés")
cluster.kmeans2 <- kmeans(mes.iris.scaled, center=2)$cluster
plot(iris.scaled.acp, col.ind=cluster.kmeans2, title="kmeans en 2 groupes")

par(mfrow=c(1,1))

Table de confusion et calcul de performances

conf.kmeans <- table(variete2, cluster.kmeans2)
kable(conf.kmeans)
1 2
notSetosa 97 3
setosa 0 50
TP <- conf.kmeans[1,1]
FP <- conf.kmeans[1,2]
FN <- conf.kmeans[2,1]
TN <- conf.kmeans[2,2]
P <- TP + FN          # nb positif dans la réalité
N <- TN + FP          # nb négatif dans la réalité
FPrate <- FP / N      # = false alarm rate
Spe <- TN / N      # = spécificité 
Sens <- recall <- TPrate <- TP / P      # = hit rate ou recall ou sensibilité ou coverage
PPV <- precision <- TP / (TP + FP)
accuracy <- (TP + TN) / (P + N)
F.measure <- 2 / (1/precision + 1/recall)
performance <- c(FPrate, TPrate, precision, recall, accuracy, F.measure, Spe, PPV)
names(performance) <- c("FPrate", "TPrate", "precision", "recall", "accuracy", "F.measure", "Spe", "PPV")
kable(performance, digits=3)
x
FPrate 0.057
TPrate 1.000
precision 0.970
recall 1.000
accuracy 0.980
F.measure 0.985
Spe 0.943
PPV 0.970
clues::adjustedRand(as.numeric(variete2), cluster.kmeans2)
     Rand        HA        MA        FM   Jaccard 
0.9605369 0.9204051 0.9208432 0.9639434 0.9302767 

Versicolor vs !Versicolor

Visualisation

variete2 <- rep("notVersicolor", 150)
variete2[variete=="versicolor"] <- "versicolor"
variete2 = factor(variete2)
table(variete2)
variete2
notVersicolor    versicolor 
          100            50 
par(mfrow=c(1,2))
plot(iris.scaled.acp, col.ind=variete2)
cluster.kmeans2 <- kmeans(mes.iris.scaled, center=2)$cluster
plot(iris.scaled.acp, col.ind=cluster.kmeans2)

par(mfrow=c(1,1))

Table de confusion et calcul de performances

conf.kmeans <- table(variete2, cluster.kmeans2)
kable(conf.kmeans)
1 2
notVersicolor 50 50
versicolor 47 3
TP <- conf.kmeans[1,1]
FP <- conf.kmeans[1,2]
FN <- conf.kmeans[2,1]
TN <- conf.kmeans[2,2]
P <- TP + FN          # nb positif dans la réalité
N <- TN + FP          # nb négatif dans la réalité
FPrate <- FP / N      # = false alarm rate
Spe <- TN / N      # = spécificité 
Sens <- recall <- TPrate <- TP / P      # = hit rate ou recall ou sensibilité ou coverage
PPV <- precision <- TP / (TP + FP)
accuracy <- (TP + TN) / (P + N)
F.measure <- 2 / (1/precision + 1/recall)
performance <- c(FPrate, TPrate, precision, recall, accuracy, F.measure, Spe, PPV)
names(performance) <- c("FPrate", "TPrate", "precision", "recall", "accuracy", "F.measure", "Spe", "PPV")
kable(performance, digits=3)
x
FPrate 0.943
TPrate 0.515
precision 0.500
recall 0.515
accuracy 0.353
F.measure 0.508
Spe 0.057
PPV 0.500
clues::adjustedRand(as.numeric(variete2), cluster.kmeans2)
      Rand         HA         MA         FM    Jaccard 
0.53995526 0.07211421 0.07722223 0.57895580 0.40737752 

Contact: anne.badel@univ-paris-diderot.fr